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COASTAL AND INLET PROCESSES NUMERICAL MODELING SYSTEM FOR 
OREGON INLET, NORTH CAROLINA 


PART I: INTRODUCTION 


Background 


1. Oregon Inlet is a large inlet through the Outer Banks barrier island 
system of North Carolina. It is located about 85 miles* south of Cape Henry, 
Virginia, and about 40 miles north of Cape Hatteras in Dare County, North 
Carolina (Figure 1). The inlet is the only existing breach in the barrier 
island system of North Carolina between the Virginia—North Carolina state line 
and Cape Hatteras. The shorelines immediately north of Oregon Inlet comprise 
a portion of the Cape Hatteras National Seashore Recreation Areas, and the 
shorelines immediately south are administered by the US Fish and Wildlife 
Service as the Pea Island Wildlife Refuge. 

2. The Manteo (Shallowbag) Bay, North Carolina, project was authorized 
by Congress in 1970 and included provisions for the stabilization of Oregon 
Inlet with dual jetties and the deepening of the ocean bar channel to 20 ft. 
In addition, there were provisions for bypassing across the inlet sand inter- 


cepted by the jetties. 


Purpose 


3. The purpose of this study was to develop a set of numerical models 
called Coastal and Inlet Processes (CIP) Numerical Modeling System to simulate 
coastal and inlet processes in the vicinity of Oregon Inlet, North Carolina. 
The models would be required to handle wave, current, and sediment transport 
processes of importance in the area and be computationally efficient enough to 


allow simulations of practical engineering utility. 


* A table of factors for converting non-SI to SI (metric) units is presented 
on page 5. 


PART II: NUMERICAL MODELS 


Introduction 


4, Although in recent years there has been considerable interest in the 
use of numerical models to simulate coastal processes, models have not been 
developed that can handle all of the complexities of an actual inlet system. 
In addition, existing numerical models that consider parts of the overall 
problem can consider only small idealized problems and not actual spatially 
large and complex problems. 

5. In order to simulate coastal processes, models must be able to con-— 
sider the propagation of waves over a complex bathymetry, the generation of 
wave-induced currents (littoral and rip currents), the circulation of tidal 
and storm surge currents, the littoral transport of sediment, and the response 
of profiles to wave action by the onshore-offshore transport of sediment. The 
models must be able to simulate processes over a relatively large area and 
resolve strong gradients in certain areas (e.g. within surf zone areas). 

6. The models described in this report use the finite difference method 
for computations. In order to cover a large region but still maintain high 
resolution in desired areas, the models use a smoothly varying grid that al- 
lows cells to be small in certain areas (e.g., surf zone or inlet) and large 
in others (e.g., ocean or sound). A piecewise reversible transformation 
(analogous to that used by Wanstrath 1977) is used independently in the x and 
y directions to map the variable grid into a uniform grid used in the computa-— 


tional space. The transformation has the following form: 


c 

me 2 + ba” (1) 
c 

y=a_+b ae (2) 


where a Des eo aq" a0 and ea are arbitrary constants for regions p 


* For convenience, symbols and abbreviations are listed in the Notation 
(Appendix A). 


and q 


nates in the computational space. 


to be centered in the computational space. 


in the x and y directions, respectively, and Oy and 5 


are coordi- 


This transformation allows all derivatives 


Many stability problems commonly 


occurring in variable grid schemes are eliminated when using this transforma- 


tion since the real space grid is smoothly varying with the coordinate and its 


first derivative being continuous at the boundaries between regions. 


Vo 


integration on a grid of spatial points. 


used with the x-coordinate increasing 


The partial differential equations are solved by finite difference 


A right-handed coordinate system is 


in the offshore direction and the 


y-coordinate increasing along the shoreline with the ocean to the right. The 
partial derivative of an arbitrary variable s in domain p is 
c) 1 
Sxeeeubetas (3) 
x i 
where 
dx ec -l 
n= 3a, = acne (4) 
Similarly 
) it @ 
ae aa (5) 
of y 92 
where 
ay cea 
Be = 2a, eae (6) 


If the grid in the x-, y-coordinate system is to have 


values of uy and UL 


constants a,b,c ,a,b_, and 
P P Pp q 


q 


of Hy and My 


will be constant (1 if 


even grid spacing, all 
Ao, = Ax and Aa, = Ay). The 


for all the domains and the values 


at grid cell faces and centers are determined using an 


interactive computer program called MAPIT. 


8. Figure 2 shows the variable 
all of the numerical models presented 
inlet processes at Oregon Inlet. The 
of approximately 60 square miles with 
as 100 ft. 


170,000 grid cells would be required. 


finite difference numerical grid used by 
in this report to calculate coastal and 
grid has 4,158 cells and covers an area 


grid cells having side lengths as small 


If a uniform grid were used with 100-ft grid cells, approximately 


Since the computational time require- 


ments of the numerical models presented in this report generally increase with 


HT il 
il 
a | 


i 
AN i 
itil 


Hit 


HT 
Mil | TMtcae oe 
aA NT 


Figure 2. 


=== 
== — 
eS) ==] Se ee Le) 
Ss en ESS Sa) 
[aa] aera] aa | a (are 
SS a baa ee er De 
= ee) a eee eS] 
(ee ie ls re ed ee) 
ES a en ea ee ee 
[= eal (eae a a | SCALE 
Hp 
a aa a a 0 5,000 10,000 FT 
sean i 
ee es 
a eas er) 
ical ef. Seed 
a eae a 


Numerical grid for Oregon Inlet simulation 


10 


the number of grid cells to approximately the 1.5 power, a regular grid with 
100-ft cells would require almost 250 times as much computational time as the 
variable grid shown in Figure 2. 

9. The grid shown in Figure 2 covers a portion of Pamlico Sound, approx- 
imately 12 miles of ocean coastlines on Bodie and Pea Islands north and south 
of Oregon Inlet, and a section of the Atlantic Ocean extending out to the 
60-ft mean low water (MLW) contour. The ocean boundary of the grid is far 
enough from the inlet to be in a region where the contours are straight and 
parallel. Thus all the effects of bathymetric complexities on waves and cur- 


rents are included in the grid computations. 


Wave Propagation Model 


10. The wave climate at Oregon Inlet is available through the Waterways 
Experiment Station Wave Information Study (WESWIS). Twenty-year hindcast data 
were obtained at a water depth of 60 ft off the coast of Oregon Inlet. The 
development of a wave propagation numerical model was necessary to propagate 
waves over the complex bathymetry around the inlet and to determine wave- 
current interactions in the vicinity of the inlet. 

11. The traditional method of propagating waves over the nearshore re- 
gion is wave ray tracing, in which the paths of individual wave rays are pro- 
pagated from deeper water to shore. One disadvantage of ray tracing is the 
difficulty of converting a large number of wave rays into arrays of wave 
heights, wave numbers, and directions of propagation for specific grid points. 
Such arrays are needed by the numerical models that calculate wave—-induced 
currents and sediment transport. It is difficult to convert a wave ray field 
into arrays providing information at specific grid points because it is not 
possible to know a priori where individual rays will propagate and thus where 
rays must be started in order to go through given points. Another disadvan- 
tage of wave ray tracing is that wave-current interaction calculations cannot 
be made, since ray theory assumes that energy cannot cross wave orthogonals. 
This eliminates the possibility of currents sweeping wave energy across 
orthogonals. 

12. In order to eliminate the problems of wave ray tracing, a numerical 
model was developed to calculate wave propagation on a gridded system. The 


model considers the same system of differential equations considered by Noda 


11 


et al. (1974) and Ebersole and Dalrymple (1980). However, to improve stabil- 
ity, some of the solution techniques employed to solve these differential 
equations were different from those used by previous investigators. The 
models used by Noda et al. (1974) and Ebersole and Dalrymple (1980) typically 
had stability problems when the angle of incidence of the waves was large or 
the bathymetry was complex. The model described in this report was developed 
also to perform calculations on a uniformly variable grid. This allowed the 
wave calculations to be performed on the same grid as wave-induced current and 
sediment transport calculations. 

13. In order to determine the angle of wave propagation 9 , use is made 


x 
of the irrotationality of the wave number vector k 


dk Os 
——— i () 
By ax (8) 
where 
k. = k cos 8 (9) 
x 
k = k sin 0 (10) 
of 
k= |k (11) 
and 09 is the angle of wave propagation. Equation 8 becomes 
08 ; 08 1 ok ; IL Os _ 
cos 8 > + sin Q me SOS az By fsin Ove 0 (12) 
14. In order to conserve wave frequency when a current is present 
w =uot+kov (13) 


fo) 


where W is a constant equal to 20/T, ; ae is the period of wave when 
there is no current, w is the radian frequency as it appears to a stationary 


observer, and ¥ is the current velocity vector. Substituting for w , using 


12 


Equations 9 and 10, and defining ¥ = Ui+ Vj yields 


1/2 


v= [ex tanh cxa)| + Uk cos 8 + Vk sin 0 (14) 


where g is the acceleration due to gravity, and d is the local water 
depth. 
15. Taking the differential of Equation 14 yields 


1 ok _ 
k Ox 
My cos 0 + tL sin 0 
62 QW gelin OS Wi cos @))_ Ox ox (15) 
ox A A 
fi Cw, - Uk cos 6 —- Vk sin 2) aa 
A sinh (2kd) ox 
where 
1 2kd % 
a3 + saan |(2- ¥ co @ — V sin s)+v cos 6+ V sin 0 
.Likewise 
—cos 0+ ol sin 0 
plipokys 00) (UNsiny Olav icosnd) Q dy 
k dy oy A A 
(w - Uk cos 06 - Vk sin 86) 
hy |e Ho) oes Oe fo. EC (16) 
A sinh (2kd) oy 
16. By substitution of Equations 15 and 16 into Equation 12 
36 sin 6 (U sin 6 — V cos 6) 08 
—| cos 0 ltt 
ox A oy 
E cos 8 (U sin 8 - V cos | 
xlsin 0 = 
A 
(17) 
ay cos 0 + au sin 0 w —- Uk cos 8 — Vk sin 0 
+ cos 8 BAG Si iS LOY BORD NO an gets Se eee ele Ue OO) 
A A sinh (2kd) oy 
cu) cos 8 + My sin 0 w - Uk cos 8 — Vk sin 0 
- sin 0 Ea ad) Sea + ee OE eT ee ese = 0 
A A sinh (2kd) Ox 
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Height equation 


17. The conservation of energy is given by 


ov. 


a || aes, eee = 
dt ox i feat aL Aa Ox 


where 


E = wave energy density and is equal to og” /8 


p = density of water 
H = wave height 


S: . = "radiation stress" 
aL) 


as defined by Longuet-Higgins and Stewart (1964). 


terms of H and assuming a steady-state wave field (dE/dt 


2 6 cos th) 6 “8 2 oan) 2S Wee 
ox g dy g 2 ax 
H 9 H Ww = w= 
YD Oy (V + C_ sin @) 3,2 yx Oy xy ox 
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S 
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S57 oy E Ez > sin 28 
where 
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Finite difference equations 


Writing Equation 18 in 


| 


(18) 


(19) 


(20) 


(21) 


(22) 


(23) 


18. Equations 3 and 5 are substituted into Equation 17, and all terms 


are written in centered finite difference form. 
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Ones read cin © ( si 0 = WY eos OD On gen Pa gen 
cos 8 L——————————— SSS 
Dy Ka rn 


cos 6 (U sin 6 —- V cos 86) 


xi sin 8 - + cos 6. 
A ; ’ 
1,J 
(24) 
a cos 0 + oN sin 0 w -—- Uk cos 8 —- Vk sin ) 
x oy ov + 2 ed - sin 06 
A sinh (2kd) oyj. . i,j 
1,J 
soll cos 0 + ou sin 7) w —- Uk cos 6 — Vk sin 8) 
ss ox ox i fo) od 0 
A sinh (2kd) ox 
i,j 
where Aa, and Aa, are the constant grid spacings in A) » G) space. 
All values of d, U , and V are known at the start of computation for 


8. If all values of ® on rows i and i+l are known, then Equation 24 
can be solved for all values on row i-l . When this has been done, the next 


row can be solved. To find 6. . it is first necessary to find k, . 
i-l,j 15J 


Since 8. j is already known, Equation 12 can be solved by Newton-Raphson 
9 


iteration for k 


53) 7 
19, Substitution of Equations 3 and 5 into Equation 19 yields 
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2g eh 6 ab 55] Mah od oo & ab 6 5) 
r) 
n | 2 (U + C_ cos 0) 5 (V + C_ sin 6) (25) 
2 dx g 2 oy g 
H G au — ow e-—- owe — ov 
+>=\o >~— +0 = = 0 
2 XX OX yx oy xy ox yy dy, Ans 


At this point in the computation, all values of U,V, d and @ are 
known, making possible the calculation of all values of k,n, g and 
Cy as well. If all of the values of H in rows i and i+l are known, 
then all values of H in row i-l can be found, which then allows the 


determination of values of H in the next row. The values of H in the two 
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outermost rows of the model grid are specified initially from wave input data, 
allowing computation of values on all other rows. 
Model use 

20. Wave climate data were provided at the boundary of the computational 
grid (Figure 2) by WESWIS. With the wave height, period, and direction of 
propagation prescribed at the ocean boundary, the wave propagation model 
determines wave heights, directions of propagation, and wave numbers at all 
4,158 grid points in the computational grid. If the wave height at any grid 
cell is greater than 0.78 of the local water depth, the model assumes that the 
wave is breaking and, based upon nonsaturated wave—breaking theory, sets the 
breaking wave height equal to 0.78 of the water depth. 

21. Figures 3, 4, and 5 demonstrate the application of the wave model to 
Oregon Inlet (existing conditions). Each of these figures corresponds to a 
particular wave condition in deep water and shows the region near the inlet 
channel and shoals. The wave height and the period in deep water are the same 
for the three cases and are equal to 2.175 ft and 6.0 sec, respectively. The 
wave direction in deep water is given by y = 45, 0, and -45 deg, respec-— 
tively, for the three cases where on is the angle the wave propagation 
direction makes with the normal to the shoreline. A uniform grid with Ax = 
Ay = 240 ft was used for these runs. In each figure, the depth contours are 
represented by dashed lines, the wave height contours by solid lines, and the 
wave directions by vectors. Figure 3 represents waves propagating approxi- 
mately in the direction of the channel. The waves are refracted onto the 
shoals, on either side of the channel, where they eventually break. As a 
result, very little wave energy is propagated up the channel. Figure 4 shows 
the inlet response for waves incident normal to the shoreline. Again the 
waves refract onto both shoals and break. There is a decrease in wave height 
just outside the shoals. Figure 5 illustrates a condition where the waves are 
obliquely incident, in a direction approximately normal to the channel. Waves 
converge on the near shoal, but note how they refract around the far shoal and 
eventually converge with reformed waves propagating across the channel. The 
computational requirements for the wave model were very modest. For the 


4,158-cell variable grid, the cost per run was under $4.00. 
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Wave-Induced Current and Setup Model 


22. When waves break, they generate currents (e.g., littoral and rip 
currents) and changes in the mean water level (setup and setdown). Since 
these currents are the main transport mechanism for sediment on a coastline, 
they must be simulated in detail in order to model coastal processes. The 
theory of the generation of wave-induced currents was developed by Bowen 
(1969), Thornton (1970), and Longuet-Higgins (1970). Numerical models have 
been developed to determine wave-induced currents (Noda 1974, Birkemeier and 
Dalrymple 1975, Liu and Lennon 1978). These models typically either consider 
only simple and idealized situations, such as plane beaches and periodic 
bathymetries, or neglect terms of the governing equations involving unstead- 
iness, advection, and/or lateral mixing. 

23. In recent years, Ebersole and Dalrymple (1980) developed a wave- 
induced current model that solves equations that include terms for unsteady 
flow, advection, and lateral mixing. The model was applied to fairly small 
problems with relatively simple bathymetries. It used a simple explicit fi- 
nite difference computational scheme and grid cells of uniform size. Sta- 
bility was obtained by Ebersole and Dalrymple by using a time-step such that 
the Courant number was less than 1.0. 

24. Since the Oregon Inlet region that required modeling was relatively 
large, it was important to develop a model that had a variable grid and was 
extremely efficient computationally. One solution technique that is extremely 
efficient is the alternating direction implicit (ADI) finite-difference 
method. ADI schemes are not limited (as are explicit schemes) to a Courant 
number less than 1 to maintain stability. Courant numbers of 5 to 10 or 
higher are typically used. In view of the similarity between the equations 
that govern wave-induced currents and currents produced by long waves (e.g., 
tides), a wave-induced current model was developed in this study by modifying 
an existing, well tested Waterways Experiment Station (WES) long-wave 
numerical model known as WIFM (WES Implicit Flooding Model) (Butler 1980). 
WIFM is a finite-difference numerical model that employs an ADI computational 
scheme and in addition uses grid cells of variable sizes. 

25. WIFM was modified to calculate wave-induced currents and setup by 
adding radiation stress terms that are the driving mechanism for wave-induced 


currents. In addition, the friction and mixing terms used in WIFM were 
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modified to conform to the formulations normally used in wave-induced current 
models. A velocity, as opposed to discharge, version of WIFM that included 
nonlinear advective terms was used. 
Equations of motion 

26. The hydrodynamic equations used in the model for wave-induced cur- 
rents and setup may be derived from the Navier-Stokes equations (Phillips 
1969). It is assumed in the derivation that the fluid is homogeneous and 
incompressible, and the vertical accelerations are negligible so that the 
pressure distribution is hydrostatic. By integrating the three-dimensional 
form of the equations in the vertical direction and applying appropriate 
boundary conditions, the depth-averaged two-dimensional form of the equations 
of motion and continuity are obtained. These equations are derived by time- 


averaging over the wavé period. The momentum equations (Figure 6) are 


aU aU aU an 1 1 XX xy EXSY, 
gee ke ed Ae ee Bu EI a Se aH EV) 
se Uae YY ap 8 oe ae Ue” pa ee Toe 0 oy t20) 
4 
— aS aS oT 
ov OV ov ony ele ES D9 2 PRI AA) FU SEE ENES Za 
og Oe Op’ Boe’ ae So od me Yr prox 7 Ku 
The continuity equation is 
an, 2 3 i 
yet oe WU) © es (Vd) = 0 (28) 


Here U and V are the depth-averaged horizontal velocity components at time 
t in the x and y directions, respectively; n is the mean free surface dis- 
placement; p is the mass density of sea water; d = n+ h is the total 
water depth where h is the local still-water depth; Thx and Sp are the 
bottom friction stresses in the x and y directions, respectively; Shose 9 Sy? 
and Syy are the radiation stresses which arise because of the excess momen- 
tum flux due to waves; and oy is the lateral shear stress due to turbulent 
mixing. Note that when the still-water level is zero, the condition n > 0 

is called "setup" and <0 is called "setdown." 


27. The numerical model uses a linear formulation for friction (Longuet- 


Higgins 1970). Thus, 


aU (29) 
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Figure 6. Definition sketch for the wave-induced current model 
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where c is a drag coefficient (of the order of 0.01) and <u > is the 


orb 
time average, over one wave period, of the absolute value of the wave orbital 


velocity at the bottom. Equations 29 and 30 are based on the assumption that 
the velocity components U and V _ of the current are small compared with the 


wave orbital velocity. From linear wave theory, 


22, 
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loop 


where H is the local wave height, T is the wave period, and k is the 
local wave number. The numerical model can be adapted easily in the future to 
other formulations for friction such as nonlinear friction. 

28. As mentioned previously, the radiation stresses are of major impor- 
tance since they furnish the "driving" forces for wave-induced currents and 
setup. For monochromatic waves, they are defined by Equations 20, 21, and 22 
in terms of the local wave climate variables H , k , and 9 which are 
obtained from the wave propagation model described previously. 

29. In the numerical model, the coordinate scheme is chosen such that 
x ais positive in the offshore direction and y is approximately in the 
alongshore direction (Figure 6). An eddy viscosity formulation is chosen for 
the lateral shear. The eddy viscosity is assumed to be anisotropic. Denoting 
S and ey as the eddy viscosities in x and y directions, respectively, in 


general, Sy is assumed to be a constant and e, a function of x and 


y .- Accordingly, 


dU oV 
ae = o(s, By tae 4 (32) 
For field applications, the eddy viscosity oe is chosen according to the 


following relation given by Jonsson et al. (1974): 


Connie (33) 


H- pT 
anh 

This represents twice the value used by Thornton (1970). It was believed that 
for field situations Equation 33 represented the eddy viscosities more realis- 
tically than the relation suggested by Longuet-Higgins (1970) for plane 
beaches. The value of ey was in general taken to be equal to the value of 
é, at the deepest part (usually near the offshore boundary) of the numerical 
grid. The numerical model is flexible enough to permit other formulations for 
eddy viscosity in the future, as our understanding improves. 

30. Using the variable grid formulation discussed previously, the 


momentum equations become: 
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The continuity equation becomes: 
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where the subscripts t , 4% , and a, indicate partial derivatives with 


respect to time, A,» and a, 9 iin Deere and the grid expansion coef- 
ficients uy and are defined by Equations 4 and 6. Note that in 
obtaining Equations 34, 35, and 36, the assumptions made in paragraph 29 were 
used. 

31. The nonlinear advective terms in the equations of motion often pose 
stability problems. These terms are handled in the present model by using a 
special scheme which will be described in the next paragraph. The eddy 
viscosity terms can cause difficulties also during the numerical computation. 
The finite difference schemes selected in the model and the formulation for 
eddy viscosity adopted in the model minimize such difficulties and stability 
problems provided that time and space steps and eddy viscosity coefficients 
are properly selected for the phenomena being simulated. 

Computational techniques 

32. In order to solve the problem under consideration on a digital com-— 
puter, the differential equations (Equations 34-36) have to be expressed in a 
finite difference form. In the present case, an alternating direction, 


implicit, finite difference scheme is employed. In view of the presence of 
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the nonlinear advective terms, a particular implicit scheme known as the 
Stablizing Correction (SC) scheme is used. The basic idea of this scheme is 
as follows. The time level is indicated by a superscript r . The scheme 
involves variables at three time levels. The values of the variables at time 
levels r-l and fr are known from previous computations or prescribed 
initial conditions. To advance the solution from time level r_ to the new 
time level r+l , an intermediate time level solution denoted by the super- 
script * is introduced. Equations 34-36 are operated in a two-step proce- 
dure. In the first step, the rectangular grid is swept in the x(a, ) direc- 
tion, advancing the solution from time level r to * . Next, the grid is 
swept in the y(a, ) direction, advancing the solution from time level * to 
r+l . The two sweeps together constitute a full time-step At . 

33. Before details of the double sweep technique are discussed, the 
notation used for individual cells of the rectangular grid will be defined. 
Let Ax and Ay denote the cell dimensions in real space in the x and y 
directions, respectively. These dimensions may vary from cell to cell. Let 
the corresponding dimensions in computational space be Aa, and Aa, 0 
These dimensions are the same for all the cells in the grid. Let m and 
n denote indices corresponding to the center of an arbitrary cell (Fig- 
ure 7). All the variables except the velocities U and V are defined at 
the cell centers. Velocities U and V are defined at cell faces m+(1/2) 
and n+(1/2) , respectively. In the x-sweep, the x-momentum equation is 
centered about the cell face m+(1/2) , and the continuity equation is 
centered about the center of the cell. The two equations are solved, using in 


‘4 r+1 : = 
the process the result U = U - At the end of this sweep, n_ and 


ut! are known. Next the grid is swept in the y direction. In this sweep, 


the y-momentum equation is centered about the cell face n+(1/2) and the 


continuity equation about the cell center. Upon solving the two equations, 


ad and yee for each cell are obtained. Thus the two sweeps 


together complete the solution for ae ’ ue » and wee : 


—r 
the values n 


34, Even though the SC scheme has been described so far in terms of the 
x-, y-coordinate system for convenience, in reality the technique must be 
applied to the equations of motion in the computational O, » % space. 
After the application of the technique, the following finite difference 


equations result. (Hereafter the bar over n is dropped for convenience.) 
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Figure 7. Cell notation 


For the a, (x)-sweep (taken along a grid cell column parallel to the 


a,-axis), 
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+ — 6 (tz) =0O at (n,m) (38) 
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In the above equations, a single bar represents a two-point average and a 


double bar a four-point average. The difference operator oF, is defined as 


8 (2) es % = (39) 
i Sply2 SAD 


for any variable Z. The definition may be extended to the operators 859 


and 6 ; >: 
a.a. 
1 j 


35. Equations 37 and 38 may be rearranged so that the unknown quantities 


are to the left and the known quantities are to the right, as follows: 


r+1 
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36. Consider the set of cells for which the index n is constant and 
equal to N . Suppose at the upper boundary cell (m = M), the velocity 


Un M+1/2 is always known. Similarly, suppose at the lower boundary cell 


(m = L) the water level is always known. Then the set of equations for 


n 
N,L 
all the cells can be written in the following matrix form if the common sub- 


script N is dropped: 
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37. Since the first matrix on the left hand side of Equation 48 is tri- 
diagonal, the above matrix equation can be solved by recursion. In general, 


the recursion relations may be written as 
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38. Since in FORTRAN computer language fractional indices are not 
possible, a new integer index system is adopted in the program. Thus all the 
variables defined at the center and the faces m+l/2 and nt+l1/2 of a cell 
n,m will be designated by the integer indices N,M . The only exceptions are 
the expansion coefficients uy and oe which are defined at cell centers and 
faces. For these the following index system is adopted. For example, Hy 
at the center of cell n,m is designated by the index 2M-1 , whereas uy at 
the face m+1/2 is denoted by the index 2M , and similarly for Be 5 Using 


this new notation, the expanded form of the recursion coefficients for the 


a, (x)-sweep may be written as follows: 


N,M 
Pe 2 (52) 
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Using the same notation, the solution (Equations 49 and 50) may be written as 


a r+l 
ne ge oP tO (58) 
we GO a (59) 
N,M-1 M-1 NM * Sy-1 


For any given N , the recursion coefficients P , Q, R, and S are com 
puted, using Equations 52-57, in succession between the boundaries in the 
direction of increasing at, (x) ° The values of these coefficients at the 
boundaries depend on the types of boundary conditions encountered. Once all 
the coefficients for a given N have been determined, the values of n* 


and urttl 


for all the cells in the column are computed, using Equations 58 
and 59, in the direction of decreasing a, (x) ° By continuing to progress to 
the next higher value of N , the whole grid is swept in the ay (x) - 
direction. 

39. The development of the finite difference equations and the recursion 
relations for the a, Cy) sweep is similar to that for the a, (x) sweep. In 


this case, using the same notation as before, the recursion coefficients may 


be written as 
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The corresponding solution may be expressed as 
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Initial and boundary conditions 


40. An initial condition of rest was used in the wave-induced current 
and setup model. Thus Tike U , and V were zero at the start of the 
calculations. Radiation stress gradients were then gradually built up to 
their full values over a number of time-steps. The gradual buildup helped 
avoid transients caused if the system were "shocked" by the sudden application 
of radiation stress gradients. The solution was stopped when a steady state 
was reached. Although the model has the capability to allow flooding of dry 
land, wave setup was not great enough to cause flooding of the grid cells used 
for Oregon Inlet. Thus along the shoreline a "no flow" (wall) condition was 
assumed at the still-water line, and no flooding of the beach was allowed. 

For the lateral boundaries, a flux type boundary condition was used. That is, 
the flux at a boundary cell was made equal to that at the next interior cell. 

41. On the offshore boundary, it is common practice to use conditions of 
either "no flow" (wall) or constant elevation. However, both of these 
conditions are highly reflective; and, as a result, the transients that 
develop during the start-up of the numerical solution tend to bounce back and 
forth between the offshore and nearshore boundaries and take a long time to 
damp out. Figure 8 shows a typical transient developed for such a case. 

These transients are evident in the results (Figure 9) presented by Ebersole 
(1980). They produce significant problems in sediment transport models. For 
example, for fixed wave conditions a steady state should develop such that 
there is no on- or offshore movement of water along a straight coastline. If 
transients are present that reflect between the coastline and the offshore 
boundary, a steady condition of no flow in the on- or offshore direction is 
not achieved. If steady state is assumed, cells will have small but steady 
currents in on- and offshore directions that will produce steady erosion and 
deposition. 

42. To eliminate the problem with transients, a radiation boundary con- 
dition was used in the wave-induced current and setup model. This condition 
allowed transients to propagate out of the computational grid. 

43. The radiation boundary condition was tested by considering a simple 
plane beach. Ebersole (1980) noticed oscillations in the wave setup at the 
beach (Figure 9), in the velocity in the offshore direction (Figure 10) and in 
the longshore velocity (Figure 11). Considering similar wave conditions, the 


model described in this report was run with the radiation boundary condition. 
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Figure 8. Setup and free gravity wave 
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Figure 9, Wave setup (Ebersole and Dalrymple 1980) 
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Figure 11. Longshore velocity (Ebersole and Dalrymple 1980) 
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Wave setup reached a steady-state (non-oscillating) value after approximately 
7-8 min (Figure 12), the velocity in the offshore direction reached a steady 
zero value (velocity should be exactly zero) after approximately 7-8 min 
(Figure 13), and the longshore current reached a steady value after 11-12 min 
(Figure 14). These tests demonstrated how well the radiation boundary condi- 
tion worked. 

Model tests 

44, The wave-induced current and setup model was tested by considering 
plane beach problems where either laboratory data or analytic solutions were 
available. The first test considered wave setup and was run for a case of 
normal wave incidence on a plane smooth laboratory beach. The laboratory 
tests were conducted by Bowen, et al. (1968). Waves with a period of 1.14 sec 
and deepwater wave height of 6.45 cm were propagated on a 1 on 12 beach slope. 
The numerical model used a 50 by 3 variable rectangular grid which represented 
the entire 40-m-long laboratory flume. Walls were used for the lateral 
boundaries as well as the offshore boundary in order to correspond to the 
laboratory situation. The solution allowed for the effect of setup on the 
wave heights in the surf zone. The solution was an iterative process since 
the setup depends upon wave heights which in turn depend upon setup. As the 
solution proceeded, the wave heights for cells in the surf zone were computed 
for each time-step by using H = y(h + n) , where n is the setup and y is 
a breaking index. The radiation stresses were changed accordingly. As sug- 
gested by Bowen, et al. (1968), a y of 1.15 was used. A spin-up time of 
10At was used. 

45. A comparison of the steady-state setup values from the numerical 
model (after 150 At) with those observed by Bowen, et al. (1968) is shown 
in Figure 15. There is excellent agreement in the offshore region. In the 
surf zone, the numerical model predicts higher setups than observed. This is 
not surprising since the numerical model did not allow flooding and runup. 

The slope of the mean water line in the surf zone is very similar in the 
laboratory measurements and numerical model results (Figure 15). 

46. A second test case considered was that of waves approaching a plane 
beach at an angle. The deepwater angle of approach 9 was 20 deg, the deep- 
water wave height was 10 ft, and the period was 12 sec. A plane beach with a 
constant bottom slope of 1 on 30 was used. A drag coefficient c of 0.01 and 


a breaking index y of 0.82 were used. Lateral mixing and advection were 
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Figure 15. Comparison of numerical solution for 
setup with experimental data 


neglected in the calculations to allow comparisons with an analytical solution 
of Longuet-Higgins (1970). Using a 100 by 6 uniform grid with Ax = Ay = 
60 ft, and a time-step of 0.5 sec, a steady-state longshore velocity distri- 


bution was obtained after a time of 800 At . Next, the same case was rerun 
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with a finer grid (Ax = Ay = 30 ft) and a time-step of 0.25 sec. Figure 16 
shows a comparison between the numerical model calculations for longshore 
current and the analytical solution of Longuet-Higgins (1970). As the cell 
size and time-step are reduced, the numerical solution approaches the 
analytical solution of Longuet-Higgins. The grid cell size and time-step 
would have to approach zero for the solutions to agree identically, since the 
Longuet-Higgins analytical solution has an infinite gradient at the breaker 
line that only can be resolved with infinitely small grid cells. Of course, 
when lateral mixing is considered, the infinite gradient disappears. Figure 
17 shows the effect of lateral mixing on the numerical solution for longshore 
velocity distribution. The mixing parameter P was defined by Longuet- 


Higgins (1970) as 


oe (68) 
where tan 8 is the bottom slope of the plane beach and N is an empirical 
coefficient which varies between 0 and 0.016. P =O corresponds to no 
lateral mixing. Figure 17 presents the numerical solution for P between 
0.01 and 0.4. For completeness, the analytical solution of Longuet-Higgins 
for P= 0 is also shown in the figure. Based upon laboratory data, Longuet-— 
Higgins suggested that P varies generally in the range of 0.1 to 0.4. 
Increasing values of P (therefore increasing lateral mixing) reduce the 
magnitude of the peak velocity, move the location of the peak velocity closer 
to the shoreline, and increase the velocities offshore of the breaker line. 
Figure 17 demonstrates that the numerical solution does indeed exhibit the 
proper behavior. 

47. To demonstrate an application of the wave-induced current and setup 
model to Oregon Inlet, a particular wave condition that occurred during the 
1962 Ash Wednesday storm was selected. From the WESWIS, a wave with a height 
of 11.39 ft, period of 8.0 sec, and an angle 06 of 51.1 deg in 60-ft depth 
was selected. This particular wave condition was selected since the angle of 
incidence is large, and thus it is a good test of the stability of the model 
for large angles of incidence and complex bathymetry. The wave propagation 
numerical model was used to calculate the wave height, direction of propaga-— 


tion, and wave number at every grid cell. 
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Figure 17. Effect of mixing parameter P on the numerical 
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48. A time-step of 18.0 sec and a drag coefficient c of 0.01 were used 
in the simulation. Figure 18 shows the bathymetry at Oregon Inlet. A steady- 
state condition was reached after a time of 67At . Figure 19 shows the mean 
water levels (setup and setdown) and Figure 20 the velocities calculated by 
the model. The velocity vectors of Figure 18 are plotted just for every other 
cell in each coordinate direction to reduce the number of vectors and thus 
make velocity patterns more apparent. In addition, the plotting of velocities 
with magnitudes less than 0.1 ft/sec is suppressed. 

49. Away from the inlet shown in Figures 18 and 19, the shoreline and 
the contours are approximately straight and parallel. Thus, there is a small 
setdown in the offshore area that is followed by a large setup. This is the 
expected setup and setdown pattern for a plane beach. The velocities shown in 
Figure 20 are mainly alongshore, and the velocity distribution is similar to 
that for a plane beach. 

50. The setup, setdown, and velocity patterns are more complicated in 
the region of the inlet (the central part of the grid). Here the breaker line 
is farther offshore. The depth in the main channel decreases first and 
increases later in the direction of the inlet. Because of these factors, the 
water sets up around the inlet and tends to create a flow into the inlet 
through the various channels. A part of the main alongshore flow goes around 
the channels and shoals and thus bypasses the inlet. 

51. Near the shoals, the patterns of mean water level and velocity are 
irregular. This is because the waves refract around the shoals and break, 
creating local setups and currents that do not necessarily conform to the 
general pattern. As the waves go toward the barrier islands, they sometimes 
re-form after breaking because the depth increases. 

52. The central processor unit (CPU) time for running the Oregon Inlet 
grid to a steady-state condition was approximately 15.5 sec on a CRAY-1 
computer. The total cost for the job, including program compilation, CPU 
time, and data file manipulation, was approximately $10. Thus the wave- 
induced current and setup model is sufficiently efficient that it can be 
applied to large practical coastal and inlet problems. If instead of an 
alternating-direction implicit finite difference model that uses a variable 
grid (the model described in this report) an explicit finite difference model 
that uses a uniform grid were used in computations for Oregon Inlet, the 


computational requirements would have been approximately 3,000 times greater. 


41 


Re 
Mi 


\\ 


\\W: * 


ane ey 

BREN \ 
RW NOR) 
AWAY 


Figure 18. Topography used for Oregon Inlet 
numerical model simulation 


°, 
enniteg:? <8) 
SEAN 


Figure 19. Mean water levels from Oregon Inlet 
numerical model simulation 


42 


An phs BAS 
KKETOTE SS F 
x<<<<4 YSN A 
4 ~ 


xX re “<<<<< < 
A Anweo<<<< 
— A Vz 


; vv 7 
= = V> 
a 
<<<<<<<e Sn Se eS 
SLAY y+ MY Y<SS 
a 


S 
<é. 


> 
SHAKE PY 
\ BoA 

IN 
pV DST 


Zz 
A 
L 
\ 


—— VELOCITY = 5 FT/SEC 


Figure 20. Velocity vector plot for Oregon Inlet simulation 


Tidal and Storm Surge Models 


53. Tidal and storm surge elevations and currents were determined using 
the WIFM model. These calculations were made in a separate WES study 
(Leenknecht, Earickson, and Butler 1984), and the results were provided for 
this study. 

54. In order to perform tidal and storm surge calculations, a larger 
area must be modeled than that modeled by the CIP numerical modeling system. 
Figure 21 shows the limits of computational grids used for the tidal and storm 


surge computations. In addition, the limits of the shore processes grid used 
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Numerical grids used for tidal and storm surge computations 


in this study are shown. The storm surge ("offshore model") grid covers a 


very large area since it is necessary that the winds that produce the surge be 


modeled. Wind fields associated with storms can cover very large areas. The 


tidal grid ("nearshore model") covers a smaller area than the storm surge grid 


but a larger area than the shore processes grid. The tidal grid must cover a 


moderately large area in order to include all of the inlets that connect the 
sound areas with the ocean. 


55. Tidal and storm surge computations performed on the grids shown in 


Figure 21 were used to provide boundary conditions for tidal and storm surge 
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computations performed on the shore processes grid. Figures 22 and 23 show 
typical tidal computations performed on the nearshore processes grid. The 
figures show mean ebb tidal currents for existing conditions and for two 
jetties with a 2,500-ft spacing. Only flows in the vicinity of the inlet are 
shown. Part of the ebb-tidal delta shoal was exposed at this time in the 
tidal cycle. 

56. The numerical tide and surge model for Oregon Inlet (Leenknecht, 
Earickson, and Butler 1984) was used to simulate also storm surge elevations 
and currents for two storms of record, namely, Hurricane Donna (1960) and the 
1962 Ash Wednesday storm. Results obtained for the 1962 Ash Wednesday storm 


were used in this study, as discussed in a later section. 


Sediment Transport Models 


Transport inside the surf zone 


57. Inside the surf zone it is the wave breaking process that is pri- 
marily responsible for the transport of sediment. This process is quite 
complex and not well understood. There is even considerable disagreement on 
the primary mode (bed load or suspended load) of sediment transport in the 
surf zone (Komar 1978). Thus a model that determines transport in the surf 
zone must be empirical, to some degree, in its formulation. 

58. The surf zone transport model developed in this study is based upon 
an energetics concept developed by Bagnold (1963) who reasoned that the wave 
orbital motion provides a stress that moves sediment back and forth in an 
amount proportional to the local rate of energy dissipation. Although there 
is no net transport as a result of this motion, the sediment is in a dispersed 
and suspended state so that a steady current of arbitrary strength will trans- 
port the sediment. Thus breaking waves provide the power to support sand in a 
dispersed state (bed and suspended load), while a superimposed current 
(littoral, rip, tidal) produces net sand transport. 

59. Komar (1974) has given the following theoretical velocity distribu- 
tion for a longshore current across the surf zone on a plane beach: 

P 


V = AX + BX i see © <2 < al (69) 


45 


| 
! 
\ 
\ 
\ 


ww TT SEN 


¢ 
\ 


\ 
‘ 


VS ST VSS ss 


4 


VY VNNWNWNA 


Noe Sat OP Sete ee Nn 


OG BObbOOAAAAAEAADEREAEAAAGasAC ae 


+ S“wovvopennenesassasssssesaagg ¢ 


U PUP prrrrrmrtiiisssieseeees & 


© SoS aeegeaaessssetseccooser eo ~ 


\ g eecceecececcsetetsscecoorse =e & 


UW ETTTITE LL ee 
© 66 ddaagsatststeasegagagaaeede 9 & 
(WTI Te ee aed 
bb oboscassconessvegISPIS/Iss F 4 A 
BIOPSY A A 
OPP 
MMMM" * 
tS QLLIISILLIT_-= « 
0 ARAMA 
Mii 
a Ahi“ 
1 Si fim 


SCALE 
(FT/SEC) 
; _— 

05 


VOSSSSSS 
i 


Figure 22. 


Ue 
/ 
iy 
Yi 
7 
7 
We 
2 


. 
\] 4 Se antes ee & = 
QQQuuesamesannssene & © 


QVAGRRAR RE OO 


= Woke Se ees 
4 
‘ 
| 


SSS ae 
16 


a 


Wy 410 04EARAA SS 
PMovrcereteere*~ % 4 


SS SSN S ese f 4 


a 
a 


e00bbsssssseer™® % 
CYYYYYYYYVY YY Soka 
Adan edEEBUVIWre S 
AbdLdEDDDWWwwTSs SO 
AbdeVwewewssresss SO 
IY, Sn poeenttetid > 


ee Re TON 
oe 


a 
eo 


A DPPODPP, 
aor PPP oP 


46 
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Ebb currents for 2,500-ft jetty spacing 


Figure 23. 


where 


Wet nel oe == 
Vv x 


fo} b 


(70) 


v being the longshore current velocity at a distance x from the shoreline, 


X, being the width of surf zone, and v, being given by 


#5 tan 
oss yo vl + ¥ Veh, = sin @ cos 6, (71) 


3 9 1 
Py = Sy laa sae TP (72) 
ba 3 9 1 
Po =- Te — 16 + TP (73) 
The coefficients in the solution are 
rs, 1 
SS sas DE (74) 
(P, -1)A 
BL = ———————_ (75) 
1 PL - Po 
———— (76) 


(1 + 0.375 ) 


where tan 8 is the beach slope, P is the mixing parameter of Longuet-— 
Higgins given by Equation 68, and y is the breaker index. Using the above 
velocity distribution and noting that the total sediment transport rate Q, 
across the surf zone is given by 

x 


b 
Q, = ¢, aid (77) 
(0) 


in which C, is the average sediment concentration (dimensionless) across the 
surf zone, and d is the water depth, one may obtain the following result 


after integration: 


B 
a 2, yt. 72 tan B 7 A 1 
Q = Cx 6 YZ V1 + Y Veh, (#24) sin oF cos 8, tan 8 ( + P, rn ) (78) 
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Equation 78 has the same functional form as Equation 4-50b of the Shore 
Protection Manual (SPM) (1984) if the group velocity at breaking is 


approximated by solitary wave theory and the units are converted as follows: 


0.00642 2 L/D ip 
Q ate Iga pgH, gh, (1 + y) sin 8, cos 8. (79) 


The sediment concentration C can be obtained by equating Equations 78 and 


s 
On at) Si, 


‘ela 0.0008 pg md) (80) 


. mae Wins Bs 
E\G gree pane? 


In computing C, , the following relation for (tan B/c) is used: 


tan 8B 1.38 CS 2oS JP) 


c c P, -1 
1 1 2 1/2} 
Xp (elo cane 


melt 


(81) 


Since for a particular wave condition H, and x, are known from the wave 


model, C, is known for the surf zone. 


Transport beyond the surf zone 


60. Beyond the surf zone, waves are not breaking. Currents (tidal, 
littoral, rip) still transport sediment, but the sediment load is much smaller 
than the load in the surf zone. Waves still assist in providing power to sup- 
port sand in a dispersed state. However, there is little turbulent energy 
dissipation, and frictional energy dissipated on the bottom represents most of 
the energy dissipation. Bed load is the primary mode of sediment transport 
beyond the surf zone according to Thornton (1972). 

61. Since beyond the surf zone it is the tractive forces of currents 
(including wave orbital velocity currents) that produce sediment movement, a 
sediment transport by currents approach is taken. Again, since the complete 
physics of the problem is not completely understood, a semiempirical approach 
must be taken. In this study, the approach of Ackers and White (1973) is 
followed after appropriate modification for the influence of waves. 

62. Ackers and White (1973) studied sediment transport due to currents. 


They used the results of 925 individual sediment transport experiments to 
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establish various empirical coefficients. The approach considers both sus-— 
pended load and bed load. It is assumed that the rate of suspended load 
transport is dependent upon the total shear on the bed. Therefore, the shear 
velocity vy is the important velocity for suspended load transport. Bed 
load transport, however, is assumed to depend upon the actual shear stress on 
individual sediment grains. Ackers and White (1973) assume that this stress 
is comparable with the shear stress that would occur on a plane granular sur- 
face bed with the same mean stream velocity. Thus the mean velocity of flow 
v is the important velocity for bed load transport. 

63. Considering only currents (not waves), Ackers and White (1973) 
derived sediment transport rate in a dimensionless form. For convenience in 


practical application, this may be written as: 


Vv Vana 1G m 
Pp * Avail 
where 
S = total sediment transport rate per unit width (vertically 
integrated combined bed and suspended sediment load) (£t3/sec/£t) 
D = sediment diameter which is exceeded in size by 65 percent (by 
weight) of the total sample 
C = exp [2-86 In Y - 0.4343 (1n Y)* - 8.128] (83) 
1/3 
y =p jsis-) (84) 
2 
v 
s = mass density of sediment relative to that of the fluid 
ny = 1.0 - 0.2432 In Y (85) 
9.66 
my = RY ae + 1.34 (86) 


(87) 


(88) 


vy = kinematic viscosity of fluid 
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Pp = porosity of sediment 


For values of Y greater than 60, C, Ny » my, , and A have the values 
0.025, 0, 1.5 and 0.17, respectively. 

64. The Ackers and White formulation modified for waves has been used in 
the past to determine sediment transport within surf zone areas (e.g., Swart 
1976, van de Graaff and van Overeem 1979, Willis 1978, and Swart and Fleming 
1980). It would appear, however, to be questionable to assume that an 
approach developed to determine sediment transport by current action would be 
appropriate to handle sediment transport in an area where turbulence due to 
wave breaking is the major mechanism for placing sediment in a state that 
allows transport by currents. However, outside the surf zone where waves are 
nonbreaking, the influence of waves is to increase the velocity of flow felt 
by sediment grains. Thus it is reasonable to use in the latter region a 
current action approach that is modified to consider also the current action 
exerted on the bottom by waves. 

65. Swart (1976) assumed that suspended load was the main sediment 
transported in the surf zone. Therefore, he modified the shear velocity to 


account for waves by increasing the shear velocity as follows: 


2 2) 
1 Yo 
ea) sere and current — War ucrent aye 2D oy Te Soe)) 
where 
q 1/2 

f= C) (fw, /2g) (90) 
G 3S is les (=) (91) 

1 Ty 

1 

fw) = Jonsson's (1966) friction factor based on bed roughness TY 

us = wave orbital velocity 


Equation 89 was originally developed by Bijker (1967) and modified by Swart 
(1974a). 

66. Van de Graaff and van Overeem (1979) noted that within the surf zone 
both suspended and bed load are significant and thus concluded that Swart's 
approach was not correct. They proposed increasing the mean velocity of flow 


in addition to the shear velocity by using the following equation: 
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u 
1 fe) 
Cy yave and current ~ WM eacerte ie 2 So ivi Se) 
where 
Ew, 1/2 
/ 
E! 10h 
C, = 18 log | D (94) 
fw. = Jonsson's friction factor with D as bed roughness 


Equation 92 was developed also by Bijker (1967) and modified by Swart (1974a). 


67. Beyond the surf zone where waves are nonbreaking, the waves do not 
increase the level of turbulence since turbulence is confined to a narrow 
boundary layer by the oscillating wave orbital velocities. Since the shear 
velocity is dependent upon the intensity of turbulence and thus the total 
energy degradation rather than the net traction on individual sediment grains, 
the shear velocity is not changed by wave action. With the wave-induced tur- 
bulence confined to a narrow boundary layer and the waves propagating essen-— 
tially without energy loss, the effect of waves is to increase the traction on 
individual grains by increasing the mean velocity felt by the grains. Thus 
the mean velocity of flow must be increased according to Equation 92. How- 


ever, the shear velocity must remain uncHanged. Thus Equation 82 becomes 


(95) 


with 


(96) 
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Equations 95 and 96 were used in the present numerical study for calculating 
longshore transport beyond the surf zone. 
Erosion and deposition 

68. In this report we are considering only the transport of noncohesive 
sand sediments. These sediments have relatively large fall velocities and 
thus remain in suspension (as suspended or bed load) for only short periods of 
time unless a forcing function is continually in effect. Thus these sediments 
act differently from cohesive sediments that may move considerable distances 
before returning to the ocean floor. For cohesive sediments, the concentra-— 
tion of sediment at a particular location may depend upon more than just the 
local forcing function. It may also depend upon the concentrations at other 
locations and flow patterns. For example, large wave activity at one location 
can produce a large concentration of cohesive sediments. These sediments can 
then be carried to an area of low wave activity and still remain in suspension 
for a long period of time so that the concentration is high in an area of low 
wave activity. 

69. In the case of noncohesive sediments, if a parcel of water with a 
high concentration of noncohesive sediment (e.g., as a result of large wave 
activity) enters a region of low wave activity, the material drops rapidly out 
of suspension, and a much smaller quantity of sediment is resuspended. Thus 
in any given grid,cell (with dimensions larger than the typical horizontal 
distances that noncohesive sediments travel when they drop out of suspension), 
the quantity of sediment being transported is a function of the wave and 
current activity within the cell. If adjacent cells supply sediment to the 
cell in greater quantities than the wave and current activity in the cell can 
support, then part of the sediment deposits in the cell and only the sediment 
that can be transported by the local waves and currents within the cell are 
transported through the cell. Conversely, if the cell can transport greater 
quantities of material than is supplied by adjacent cells, then erosion must 
occur within the cell as bottom material supplies the additional required 
material. For example, if an offshore breakwater is built, sediment will 
deposit in its lee where wave activity is lower than in adjacent areas. If 
something increases wave convergence in an area, erosion will occur. 

70. Since the rate of noncohesive sediment transport within a cell is 
dependent only upon wave and current conditions within the cell, transport can 


be modeled using a "box model." In this model, the transport in and out of 
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the four faces of a cell must be considered. If more material enters a cell 
than leaves, deposition of material occurs, and if more material leaves than 


enters, erosion of material occurs within the cell. If S is the total 


out 


transport rate out of a cell and Si, is the total transport rate into a 


cell, then 
(97) 


where 
Ah = change in the bed elevation of a cell 
At = time-step 
A_ = area of cell 
Pp = porosity of sediment on bed 


71. When applied to each grid cell of the grid shown in Figure 2, 
Equation 97 determines the erosion and deposition throughout the grid. It is 
assumed that there is no net erosion or deposition within all the boundary 


eas, § and Sant can be determined once the current and wave fields 


in 
have been determined by the numerical models described earlier. Within the 
surf zone, both the sediment concentration C; and the total water discharges 


(tidal and wave-induced) are used to compute Say GG! Soyrte 


Profile Response (Onshore-Offshore) Model 


Introduction 

72. In addition to littoral and tidal transport, beach profiles respond 
to the wave climate through onshore-offshore transport of sediment. This 
onshore-offshore transport occurs at a much (orders of magnitude) slower rate 
than does littoral transport of sediment (Galvin 1983). However, although the 
gross littoral transport of sediment is very large at any point in the surf 
zone, the net deficit or surplus of sediment is usually fairly small. Con- 
versely, the gross onshore-offshore transport of sediment is small, but the 
net deficit or surplus can be equal to the gross over short periods of time. 
For example, during a storm there may be offshore movement of sediment to an 


offshore bar. Although the total amount of sediment that moves to the bar is 
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small compared to the total gross yearly littoral transport, all of the sedi- 
ment that is transported must be eroded from the nearshore area. 

73. The onshore-offshore transport of sediment is not well understood. 
Attempts have been made in recent years to establish the equations governing 
onshore-offshore transport and solve them using numerical models. However, 
none of these attempts has produced a numerical model that can be used for 
reliable quantitative predictions of onshore-offshore transport. For example, 
Wang (1981) evaluated onshore-offshore models developed in recent years by 
Sunamura (1980), Dally (1980), and Yang (1981). Wang concluded that "Beach 
profile modeling is a quite recent endeavor. It is a difficult problem be- 
cause the physical process is complicated and is not well understood. The 
three models introduced here are not at operational level and are not adequate 
for quantitative predictions." 

74. Since the equations governing onshore-offshore transport are not 
completely known, in this study a numerical model is developed that is 
strongly based upon concepts developed by Swart of the Delft Hydraulics Lab- 
oratory, the Netherlands (Swart 1974a, 1974b, and 1976). Swart's concepts 
were extended to allow the model to consider a variable datum (time-varying 
tide), a variable wave climate, and onshore transport in addition to offshore 
transport (Swain and Houston 1983, 1984a, 1984b, and Swain 1984). 

Governing equations 

75. In his conceptual model, Swart divided a normal beach profile into 
three zones (Figure 24), each with its own transport mechanism. The first 
zone is a backshore above the limit of wave runup. If windblown sediment 
transport is neglected, there is no transport in this zone. The second zone 
is a developing profile (D-profile) where a combination of bed-load and sus— 
pended load transport takes place. The dividing point between these two zones 
is the highest location that waves reach on the beach. Since the tide datum 
and wave climate vary with time, this dividing point moves with time. The 
position of maximum runup was determined empirically by Swart and is given by 


the following equation (all units are metric): 


-0.000143 Ho 188 79-93 


= = eaete seen 1s ee eee 
Bea 7650 Deo 1 - exp 50.786 (98) 
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Ds5go = median particle diameter 


H = maximum wave height in the spectrum (equal 
to twice the significant wave height) 


dt = wave period 


The third zone is the transition area seaward of the D-profile and landward of 
the point where sediment motion is initiated by wave action. Bed-load trans-— 
port is normally the only transport in this zone. The point dividing the 
lower limit of the D-profile and the upper limit of the transition area was 
determined empirically by Swart, and the depth of this point is given by the 


following equation: 


Rosny teem 
ho = 0.0063 do ©xP | 0.804 0.093 (99) 
T D9 


where do is the deepwater wave length and H, is the deepwater wave height. 


76. The basic assumption in the theory of Swart is that the D-profile 
will eventually réach a stable situation under constant wave attack. This 
stable situation implies both an equilibrium form and position of the beach 
profile. By considering many small- and full-scale tests of profile develop- 
ment under wave attack, Swart was able to develop equations that determine the 
form and position of the equilibrium profile for different incident wave 
climates. 

77. At every location "i" on the D-profile, Swart defines an onshore 
and offshore segment of the profile (Figure 24). The length of the onshore 
profile is represented schematically by the distance Lj, and the length of 
the offshore profile by Ly; . The length difference at each point i be- 
tween the onshore and offshore sections of the D-profile, (Ly - Lj); , is the 
key parameter used by Swart to characterize a profile. The value of this 
parameter when the profile is in equilibrium is defined as W; . Swart thus 
represents the equilibrium profile by a "W-curve". He defines the W,; value 


at the still-water line as W- which is given by the following equation: 
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H H 
6 fe) 3],.0.132.-0.447 fo) 
WEED Gul Fs 110) Ho 50 x 


2W 
16) (0) fo) 


46 Aa) iba. £2 tore (100) 


where 6 is the total depth of the D-profile (Figure 24). The variation of 
the ratio W;/W. over the D-profile determines the dimensionless form of the 


D-profile. Swart determined the following equation for W;/W, 


W 4 
i 7 2 1.36x10 
W = 0.7 A. tedt ae Slo S)7/ es 0)" 19) Deo A. Deo (101) 
where 
In = Ons ; 6 Aa ; 
eo eyes 2i _ the dimensionless position in the (102) 
ey ) ~ D-profile measured positively downward 


from the still—-water level 


b=1 for A >O =, i.e., below the still-water level 


1) = (0) sF@Ke A. < O , i-e., above the still-water level 


78. The transport rate (Sit (i refers to position on D-profile 


and t to time) was found by Swart (1976) to be given by the following 


equation: 
He WikSie Kane 
S =|[¢ eichewa exp [ -X,t (103) 
Vit y 554 5 
where 
; | \2-66|- 0.079 
g\D 280 eo et oadomine oe Dee oa (104) 
y tt fo) fo) 50 h 
m m 
Sa ie One ID ae Wh te Ilo 
Ghat pe iG San) =n eee 
yo 6 bm 2 1°mo 2 1° io 
(y, - y,). 
ba ot eee all ato 
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X, = (107) 
2 815 994 
Yor = I = W - Ly - be (108) 
10 10 
5. 
ba See Me (109) 
Zoek 


Subscript m refers to middepth (655 = 0.56) , and subscript o refers to 


time t=O. W, is the difference in length between the transition distance 


t 
for an equilibrium srotude and the transition distance for an initial profile 
(Figure 24). 

79. Although the Swart model is a two-dimensional onshore-offshore 
model, transport depends upon the angle at which waves attack the D-profile 
area. Oblique wave attack produces increased shear stress on the bed due to 
the presence of nearshore currents generated by the oblique wave attack. 


Swart found that (sy)4 increases when the three-dimensional effect of 


oblique wave attack is included. The increase is given by the following 


factor: 
4.5 
Me 
l+m ~ (110) 
Z eTke 
where 
my = 1.94 + 2.97 sin oF (111) 
os = angle of wave incidence at breaking 
fw 1/2 
ey = cn (76) (112) 
C, = Chezy coefficient 
fw = wave friction factor (Jonsson 1966) 
ug = wave orbital velocity at the bed 
No = 1.27 - 0.39 sin 8, (113) 


The longshore current velocity v is computed using the following expression 
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tan 8B 


= (enyae sin 8 cos 8 (114) 


80. Every time the wave climate is changed, a new equilibrium profile 
(given by the W; values) must be calculated. The existing profile at this 
particular time then becomes the initial profile used to determine the (Ly - 
L1);9 Values at all points i. The rate of sediment transport (Sy). is 
then determined at each point i and is modified to include the effects of 
oblique wave attack. This rate of transport varies with time exponentially. 
Thus, the profile begins initially to move toward an equilibrium profile at a 
fast rate, but the rate slows as the profile approaches an equilibrium 
profile. Actually, the time required to reach equilibrium is long compared to 
the time a typical wave condition exists. Thus, profiles are always moving 
toward equilibrium profiles but never fully reaching them. 

Model testing 

81. The profile response model was compared with laboratory tests of 
profile modification under wave attack and with prototype measurements of 
profile modification. Figures 25 through 28 show comparisons between beach 
profiles that were measured in experimental tests by Eagleson et al. (1963) 
and results of the profile response model. The solid lines are profiles 
calculated by the model, and the dashed lines are measured profiles. The only 
parameters varied in the numerical model calculations were those that were 
varied in the laboratory test (initial profiles, wave heights, wave periods, 
and duration of the tests). The major features of the measured profiles are 
reproduced in the numerical model results. The agreement between measured and 
calculated profiles is quite remarkable considering the great complexity of 
the fluid-sediment interaction process. 

82. Comparisons also have been made between measured profile modifica-— 
tion in the prototype and numerical model results. For example, a storm 
during the period 16-21 February 1980 at Santa Barbara, California, was docu- 
mented (Gable 1981) in detail during the Nearshore Sediment Transport Study 
(NSTS). This was a large storm that produced approximately 40 m of shoreline 
erosion. Daily profile measurements were made by the NSTS in addition to 
complete directional spectral wave data. Figure 29 shows good agreement 
between measured profiles and the profile response model simulation over the 


5-day period of the storm. The only inputs to the numerical model were the 
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MEASURED PROFILES 
(DATE SHOWN) 
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Figure 29, Comparison of calculated and measured beach 
profiles for Leadbetter Beach, California, 1980 


initial profile (16 February 1980); the sand size; hourly values of signif- 
icant wave height, period, and direction; and hourly tide levels. Tidal 
fluctuations played an important role in the mechanism of cross-shore trans— 
port (Swain and Houston 1984a). The sensitivity of the model to other 
important parameters is discussed elsewhere (Swain 1984). 

83. A comparison has been made also between measured profile modifica- 
tion during the Currituck Sand-Bypass Study (Schwartz and Musialowski 1980) 
and the profile response model calculations. This study involved placement of 
26,750 cu m of sediment on the coast near New River Inlet, North Carolina, 
using the split-hull dredge CURRITUCK. Since the dump area was relatively 


small, there were significant "end effects" so that the amount of sand in the 
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profiles was not constant. However, end effects were minimized by selecting a 
profile through the center of the dump. Wave characteristics were obtained 
during the study using Littoral Environment Observation (LEO) techniques. 
Figure 30 shows an initial profile measured after the dumping was completed 
(11 August 1976) and measured and calculated profiles at the end of the LEO 
measurements. The profile response model predicted that there would be little 
modification in the profile over this time period except for some erosion of 
the break point bar and filling of the adjacent trough. The measured profile 
confirms this numerical prediction. Figure 31 shows that the calculated and 
measured profiles differ at most by a few tenths of a metre in elevation. 

This difference is undoubtedly within the level of accuracy of the profile 
measurements. In addition, LEO measurements are not precise means of 


measuring wave conditions. 


ELEVATION, M 
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Figure 30. Comparison of calculated and measured shore-normal 
profiles for New River Inlet, North Carolina 
(dredged disposal sand movement) 
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Figure 31. Comparison of calculated and measured final 
profiles for New River Inlet, North Carolina 
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PART III: APPLICATIONS 


Ash Wednesday Storm 


84. As a verification test of the numerical modeling system's ability to 
determine sediment transport near Oregon Inlet for an extreme storm event, the 
Ash Wednesday storm of 1962 was simulated. The storm lasted from 7 to 
1l March 1962. This record northeaster caused severe beach erosion along 
Bodie Island and Pea Island on the north and south sides of the inlet, respec-— 
tively. The erosion continued until it reached high dunes. On the northern 
part of Bodie Island, the land was low and flooded. 

85. To run the simulation, the wave climate (significant wave height, 
period, and direction) in deep water was hindcast at 3-hr intervals throughout 
the storm using WESWIS. Simultaneously, the tidal and storm surge levels and 
currents during the storm were obtained using the numerical tide and surge 
model. 

86. Since the major part of the net sediment transport during a severe 
storm is in the onshore-offshore direction, only the profile response model 
was run using the information on waves, astronomical tides, and storm surge 
levels. The measured shore-normal pre-storm profiles for the two islands were 
furnished by the National Park Service (NPS), US Department of the Interior 
(DOI), and were taken a year prior to the storm. The poststorm profiles were 
taken a few weeks after the storm and were furnished by the NPS. Thus the 
measured data do not represent just the erosion due to the storm. However, 
the erosion produced by the storm was much greater than normal yearly fluctua- 
tions. The measured profiles indicate that shoreline erosion varied from 
sixty to several hundred metres with higher erosion being near the inlet. 
Figures 32 and 33 present a comparison between measured shore-normal erosion 
and numerical model results. The agreement is remarkable, considering the 


uncertainties associated with the measured data. 
Evaluation Of Nonstructural Solution 
87. As an alternative to the inlet stabilization of Oregon Inlet that 
was authorized in the overall Manteo (Shallowbag) Bay project, DOI proposed a 


dredging and disposal procedure that would not require the construction of 
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Figure 32. Comparison of calculated and measured shore-normal 
erosion for Bodie Island during 1962 Ash Wednesday storm 
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Figure 33. Comparison of calculated and measured shore-normal 
erosion for Pea Island during 1962 Ash Wednesday storm 
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jetties at Oregon Inlet. DOI administers the Federal lands adjacent to Oregon 
Inlet on which the proposed jetties must be anchored. 

88. The DOI plan attempted to minimize disturbance of park and refuge 
lands by proposing that material dredged from the inlet be placed in the near-— 
shore zone as close to the beach as possible by small split-hull hopper 
dredges. This was in lieu of hydraulically transporting the dredged material 
from the inlet along the beach via a pipeline with disposal being accomplished 
in a manner similar to that in conventional beach nourishment. The DOI pro- 
posal assumed the dredged material placed in the nearshore zone would be dis-— 
persed shoreward by wave action at a rate sufficient to prevent dredging- 
induced beach erosion. 

89. The only dredging plants that have the capability to quickly dispose 
dredged material in the shallow nearshore zone in the presence of waves are 
the CURRITUCK and ATCHAFALAYA/MERMENTAU class split-hull dredges. The split— 
hull design allows these dredges to dump their loads very rapidly and at the 
same time rise quickly to a shallow draft, thus avoiding the danger of 
striking the bottom. The CURRITUCK class dredge has a hopper capacity of 
315 cu yd, an unloaded draft of 3 ft, and a loaded draft of 8 ft. The 
ATCHAFALAYA/MERMENTAU class has a hopper capacity of 1,300 cu yd, an unloaded 
draft of 5.7 ft, and a loaded draft of 14 ft. 

90. The US Army Engineer District, Wilmington (SAW) (1983), performed a 
feasibility study for the plan proposed by the DOI. Material dredged from the 
inlet would be placed in a nearshore disposal area as shown in Figure 34. The 
disposal zone would begin at a point approximately 12,500 ft from the seaward 
end of the entrance channel. At this location, the ebb-tidal delta ceases to 
exist; and bottom contours are essentially parallel to the shore. Placement 
of material closer to the inlet than this point would cause the material to be 
in the ebb-tidal delta area. Material in the ebb-tidal delta typically would 
not be transported to adjacent beaches; thus erosion of the beaches would be 
accelerated. Much of the material also would move toward the inlet and 
quickly reenter the navigation channel. The farthest disposal point would be 
located approximately 65,000 ft from the seaward end of the ocean entrance 
channel. 

91. The split-hull dredges would place material in the nearshore zone by 
dropping it in the pattern shown in Figure 34. Dimensions of the disposal 


mounds were determined by SAW based on the length and breadth of the vessel 
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hoppers and the loaded and unloaded drafts of each vessel. Disposal depths 
were -ll ft and -17 ft MLW for the CURRITUCK and ATCHAFALAYA/MERMENTAU 
classes, respectively (based upon a 3-ft clearance of each vessel to allow for 
vessel vertical motions that would occur as a result of steep shoaling waves 
in the nearshore area). The dredges would continue the dumping pattern until 
the end of the disposal area was reached. Then the pattern would begin again. 
If the material in the first disposal mound were sufficiently dispersed to al- 
low the dredges to dump material once again, the second disposal mound would 
be located where the first mound was dumped. If not, a second row of material 
would be dumped immediately behind the first. This placement of mounds would 
continue indefinitely with attempts made to dump the material in the shal- 
lowest depth of water available (depending upon the dispersal of earlier 
mounds). 

92. A key question about the feasibility of the DOI proposed nonstruc-— 
tural solution was whether the material placed in the nearshore area would be 
dispersed shoreward by wave action at a rate sufficient to prevent dredging-— 
induced beach erosion. To address this question, only the profile response 
(onshore-offshore transport) model described in this report was used since the 
material to be deposited would be away from the area where tidal currents are 
significant and in an area where bathymetric contours are approximately 
straight and parallel. Changes in the rates of onshore-offshore sediment 
transport resulting from the disposal of the dredged material were determined 
from the differences between the sediment transport rates computed for a 
control profile (i.e. a profile without the dredged material) and the sediment 
transport rates computed with the dredged material in place on the profile. 
The calculations were accomplished using wave and tide conditions for 1975 
which were determined to represent typical conditions in the Oregon Inlet 
area. 

93. The hindcast wave conditions for 1975 were provided by the WESWIS in 
10 m of water and at 3-hr intervals throughout the year. The astronomical 
tide for 1975 was generated using tidal constituents. SAW supplied an initial 
profile that was representative of the region where the dredged material would 
be dumped (Plate 1). 

94. The profile was modified by making comparisons between the dumped— 
material profiles and the control profile every time-step. The width of the 


surf zone was determined by assuming that the breaking depth was the shoaled 
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wave height divided by 0.78. The dumped-material profile was compared with 
the control profile over this surf zone region. When the dumped-material 
profile had an excess of material relative to the control profile over any 
section of the surf zone region, the excess material was transported out of 
the dumped-material profile in the following manner. If the excess of mate- 
rial was less than the quantity of material that could be transported accord- 
ing to the longshore transport formula of the Coastal Engineering Research 
Center (CERC) given in the SPM (1984), all of the excess sediment was trans- 
ported out of the profile; however, no additional sediment transport was 
computed as would arise from erosion of surf zone and beach line. Thus, the 
focus of the computations was restricted to the fate of dredge disposed mate- 
rial on entering the surf zone and not on erosion of the existing control 
profile landward of the breaker line. Om the other hand, if the quantity of 
material entering the surf zone was greater than that which could be trans- 
ported during the time-step according to the CERC formula, then the material 
transported out of the dumped material profile was limited to the quantity 
given by the CERC formula, thereby showing surf zone accretion. It was recog- 
nized that an excess of material over just a part of the surf zone would not 
be transported by the full power available in the surf zone. Thus, it was 
assumed that the material transported in the surf zone was distributed 
linearly over the ,surf zone with zero transport at the shoreline and maximum 
transport at the breaker line. If the excess of material was distributed over 
a section of the surf zone, only that power available over the section was 
used to transport material. Tests showed that the fate of the dumped material 
did not depend significantly on the assumed transport distribution across the 
surf zone. Virtually the same results were obtained when the transport 
distribution was assumed to be uniform over the surf zone. Apparently the 
quantity of excess material in the surf zone region during any time-step was 
usually less than the quantity of material that could be transported according 
to the CERC formula so that all of the material was transported during the 
time-step. 

95. If there is to be shoreline stability, the dumped material must 
enter the nearshore region. During large storms, the surf zone can extend a 
very large distance seaward. Material in fairly deep water will be trans-— 
ported in a direction parallel to the coast during these large storms. How- 


ever, if this material in deeper water does not move toward the shore, it will 
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not aid shoreline stability since once the storm is over this material will 
still be in deeper water and not in the normal littoral regime area. To 
include this effect in the calculations, material was transported out of the 
dumped material profile over a surf zone width that covered the normal near-— 
shore transport zone. This normal nearshore transport zone was assumed to be 
contained within the surf zone which existed 95 percent of the time over the 
typical year. At Oregon Inlet, this zone extended to a water depth of 
approximately 2.6 m. 

96. There were two depths at which the dredged material was dumped: 
11 ft (3.4 m) and 17 ft (5.2 m). The shape of the dumps was trapezoidal, as 
provided by SAW (Plates 2 and 3). The dump in 11 ft of water was 11,000 ft 
long, and the dump in 17 ft was 9,800 ft long. Multiple rows of dumps were 
used. For example, an initial dump in 11 ft of water was made. The next 
month the profile was monitored to determine whether or not the material in 
the dump had moved out of the region sufficiently (no more than a 0.5-ft depth 
of material remaining) to allow another dump to be made. If a dump could not 
be made, a second row offshore from the first row was established. The next 
month the profile was monitored again. If the second row contained material 
that extended more than 0.5 ft above the base of the first dump location, then 
it was assumed that the dredge could not pass, and the material was dumped in 
a third row. The first row location was monitored if the dredge could pass 
the second row location. If the material could not be dumped in the first row 
location, it was dumped in the second row location again. A maximum of three 
tows was allowed for the dumps beginning in 17 ft of water, since it was 
desired not to dump material in water depths so great that the material would 
not be transported. Originally, the 17-f£t dump was 4,900 ft long, and mate- 
rial was dumped once a month. However, three rows of material were not suffi- 
cient for the case (near the end of the year a dredge could not pass the third 
row and thus could not dump material). Therefore, the row was extended to a 
length of 9,800 ft, and dumps were made once every 2 months (so that the same 
total quantity of material was dumped). The 11-ft dumps were made once a 
month. Several dumping schedules were tested. First, a year was simulated 
with dumps beginning in January (both 11- and 17-ft depths). Then similar 
calculations were made for years starting in April, July, and October. 
Table 1 presents the percentage of dumped material transported into the active 


surf zone in 1 year. Somewhat more material was transported into the active 
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Table 1 


Percentage of Dumped Material Transported 


into the Active Surf Zone in 1 Year 


Starting Month ll-ft Depth, percent 17-ft Depth, percent 
January 27.6 D2)Q©) 
April 2U 58 Za 
July 25.8 23.4 
October 28.8 24.8 


surf zone for the 11-ft depths than for the 17-ft depths. The quantity of 
material transported into the active surf zone did not depend significantly on 
the initial month of placement. During the tests, there was no noticeable 
tendency for material to be transported to the nearshore region and remain 
there permanently. Material was dispersed by the littoral currents when it 
entered the nearshore region. This effect was noticed during the Currituck 
experiment also. 

97. In addition to the extensive l-year simulations, an overall 5-year 
simulation was made to determine the volume rate of shoreward transport of an 
initial l-year dump of material over a consecutive 5-year time span. In this 
simulation, a l-year supply (1.45 million cu yd) of dredged material was 
placed in a nearshore water depth of 17 ft. The results of this test, in 
terms of percent of the initial volume moving into the surf zone each year, 
were as follows: (a) 1st year--25.5 percent; (b) 2nd year--27.5 percent; 

(c) 3rd year--26.5 percent; (d) 4th year--18.8 percent; and (e) 5th year--1.4 
percent. From this analysis, SAW adopted an overall annual shoreward volume 
rate of transport of 25 percent of a single year's volume of disposed material 
and evaluated the shoreline response for this rate of material transport. SAW 
found that this rate of onshore movement of dumped material was inadequate and 
would result in severe erosion occurring along the northernmost 3 miles of Pea 
Island. Accordingly, SAW recommended that no further consideration be given 
to the DOI dredging/nearshore disposal plan. 

98. SAW contracted with Coastal and Offshore Engineering Research, Inc. 
(COER) to perform the same basic analyses described above, except an independ-— 
ent numerical model (Perlin and Dean 1983) developed by COER for CERC was 


used. In the COER model the bathymetry was represented by n-contour lines, 
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each of a specified depth. A requirement of the model is that depth must 
increase in the offshore direction. Therefore, the dumps placed in the DOI 
dredging plan were represented as a flat area of the bathymetry. To simulate 
the placement of dumped material at a depth of 11 ft, material was placed 
between the 7- and 1l- ft contours. Material was placed also between the 

1l- and 14-ft contours. 

99. The COER model predicted that between 15 and 35 percent of the mate- 
rial added between the 7- and 11-ft contours, and between 5 and 25 percent of 
the material added between the 11- and 14-ft contours would be dispersed into 
the nearshore transport system during the first year. These percentages of 
material that would move ashore during the first year are similar to those 
calculated by the model described in this report. Quantitative comparisons . 
between the models are not possible since the COER model cannot simulate the 


actual bathymetry of the dumping plan. 
Wave-—Current Interaction 


100. For cases with and without the jetties, SAW performed an ocean bar 
channel dredging analysis. A needed input to the analysis was the period of 
time that split-hull hopper dredges of the CURRITUCK and ATCHAFALAYA/MERMENTAU 
classes could operate in the channel under the influence of waves. Waves tend 
to be large near the ocean bar channel as a result of shoaling and refraction 
on the ocean bar and wave-current interactions. The wave climate at Oregon 
Inlet was known offshore at a water depth of 60 ft from WESWIS. The wave 
propagation model described in Part II was used to transform the wave climate 
at a water depth of 60 ft to the ocean bar channel (including the effects of 
wave-current interaction). 

101. The wave propagation model was applied to Oregon Inlet using wave 
periods of 7, 9, and 12 sec, wave heights of 1 to 6 ft in 1-ft intervals, and 
the numerical grid shown in Figure 2. Tidal currents were obtained by using 
the tidal circulation numerical model described in Part II. Jetty spacings of 
2,500, 3,500, and 5,000 ft and existing conditions without jetties were 
considered. 

102. Figures 35-37 show some typical results produced by the wave propa-— 
gation model. For example, Figure 35 shows the effect of ebb current magni- 


tude on wave amplification for a 12-sec wave with an incident wave height of 
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Figure 35. Effect of ebb current magnitude on wave 
amplification (no jetties) 
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Figure 36. Dependence of wave amplification 
on wave period 
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Figure 37. Effect of jetty spacing on wave steepness 
(H/L = 0.057 for no jetties) 


3 ft (60-ft depth) and the condition of no jetties. The wave height increases 
with increasing ebb current, as expected. Figure 36 shows the dependence of 
wave amplification on wave period for a 2,500-ft jetty spacing. The wave 
amplification increases with increasing wave period, since the longer the wave 
period the sooner the wave interacts with the bottom and refracts and shoals. 
Figure 37 shows the effect of jetty spacing on wave steepness for a 9-sec wave 
with an incident wave height of 3 ft (60-ft depth). As expected, an increase 
in jetty spacing causes a decrease in the wave steepness since it reduces the 
ebb velocity. The decrease in wave steepness is not dramatic, however, since 
the jetty spacing does not change the ebb currents very much, as seen in 
Figures 38 and 39. 

103. The wave climate at a water depth of 60 ft is in terms of signif- 
icant wave height. Significant wave height is an average of the one-third 
highest waves of a given wave group. SAW determined that the controlling 
factor in vessel operation is the highest wave in a wave group. Thus SAW 
calculated the highest 1 percent wave for each spectrum represented by a 
significant wave height. Using the wave propagation model, these waves were 
then propagated to the ocean entrance channel. 

104, The limiting wave heights for the operation of the two classes of 


split-hull dredge were determined by SAW based upon information from the 
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Ebb currents for 3,500-ft jetty spacing 


Figure 38. 
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Ebb currents for 5,000-ft jetty spacing 


Figure 39, 


operators of the vessels. In the case of the CURRITUCK, operations ceased 
when wave heights encountered were in the range of 6 to 8 ft, whereas the 
ATCHAFALAYA/MERMENTAU operated up to about 10-ft wave heights. Using the 
results of the wave propagation model calculations, SAW determined limiting 
wave heights for dredging operations in Oregon Inlet to be deepwater signif- 
icant wave heights of 3.0 ft for the CURRITUCK and 4.0 ft for the 
ATCHAFALAYA/MERMENTAU classes. 

105. The wave height limitations for the dredging operations predicted by 
the wave propagation model were substantiated by SAW which compared the daily 
logs of the dredge HYDE (which operated at Oregon Inlet during 1960-1971) to 
wave data recorded during most of this period. The HYDE was a seagoing 
trailer-suction hopper dredge with seakeeping qualities similar to those of 
the ATCHAFALAYA/MERMENTAU class dredge. The wave gage was located approxi- 
mately 8 miles north of Oregon Inlet on a fishing pier at Nags Head. On the 
basis of the comparison, it was concluded that the HYDE ceased to operate when 


significant wave heights reached 4.0 ft at the wave gage. 


Single Jetty Evaluation 


106. The sediment transport model developed for transport within the surf 
zone and beyond was applied to Oregon Inlet for evaluating erosion and accre-— 
tion in the ocean bar entrance channel and the lateral movement of the bar 
channel for the case in which just the south jetty was in place. This single 
jetty condition simulates a construction sequence in which construction of the 
south jetty would be completed prior to the beginning of construction of the 
north jetty. 

107. The wave climate for an average year at a depth of 60 ft MLW near 
Oregon Inlet was obtained from WESWIS based on 20-year hindcast data. This 
information is given in the form of percent of occurrence of wave height, 
period, and direction combinations. For running the sediment model, the wave 
climate information was further consolidated into 35 combinations of wave 
height, period, and direction, and the percent of occurrence of each 
combination was determined. Table 2 shows the 35 combinations. The wave 
propagation model was run with the south jetty in place for these combina- 
tions, and the local wave information at each of the grid cells and the 


breaker line(s) was determined for each case. This information was used as 
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input to the wave-induced current model which computed and stored on files the 
steady-state wave-induced currents and setups for the single south jetty case. 

108. To run the sediment transport model for an average year, the year 
was divided into blocks of 4—hr duration during each of which the wave climate 
was assumed to be constant. During the running of the sediment model, random 
numbers between O and 1.0 were generated every 4 hr of prototype time using a 
random number generator subroutine. Based on the value of the random number, 
one of the 35 combinations of wave climate listed in Table 2 was picked, and 
the corresponding file from the wave-induced current model was used. For 
example, if the value of the random number generated was between O and 0.2580, 
the first combination listed in Table 2 was selected. If the random number 
was between 0.2580 and 0.3725, the second combination of Table 2 was picked, 
and so on. In this way when the full year simulation was run, the percent of 
occurrence of each combination in the simulation was equal to the percent 
determined from WESWIS, and the wave climate for the average year was faith- 
fully reproduced in the numerical model. 

109. Tidal currents and elevations for the mean tide were obtained by 
applying the numerical tide and surge model to the nearshore processes grid. 
The tidal elevations and currents were available for this study at intervals 
of 10 min for one full tidal cycle. An examination of these data revealed 
that the tidal elevations and currents did not vary significantly during 1 hr 
of record. Therefore, tidal elevations were varied hourly in the sediment 
transport model. In addition, to be consistent with the wave-induced cur- 
rents, tidal currents were varied every 4 hr in the sediment transport model. 
During the l-year simulation, the results of the numerical tide model were 
repeated at the end of each tidal cycle. 

110. For purposes of computation, the entrance channel for Oregon Inlet 
was defined on the basis of the maximum water depth. In addition to the cell 
having the maximum water depth in a given row of cells, two cells on either 
side were included in the channel, provided their depth exceeded 10 ft. Fig- 
ure 40 shows the channel so defined. The sediment model was modified to simu-— 
late removal of 40,000 cu yd of dredged material per year from the channel. 

111. The wave-induced current model provided the discharge on each cell 
face (longshore and cross-shore direction). The tidal velocities obtained 


from the numerical tide model were used to calculate the tidal discharge on 
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Table 2 
Selected Wave Characteristics from 20-Year 


Hindcast of WESWIS 


Wave Angle 
Wave Wave with the 

Serial Height Period Shoreline Percent of yi Percent of 
No. ft sec deg Occurrence Occurrence 
1 0.82 8.0 45 25.80 25.80 
2 2.46 UsS 45 11.45 S)/ 6225) 
3 4,10 Vos 45 Sr50) 40.75 
4 5.74 V5 45 iG abs) 41.90 
5 7.38 Use 45 0.50 42.40 
6 9.02 708 45 0.20 42.60 
7/ 10.66 8.5 45 0.05 42.65 
8 0.82 Vod 75 2.70 45.35 
9 2.46 Uoe) US 3.45 48.80 
10 4.10 1260) 75 2.05 50.85 
11 5.74 12.0 75 0.80 51.65 
12 Uods oe U5 0.35 52.00 
13 9.02 Vod) iS 0.35 DAo35) 
14 10.66 8.5 75 0.15 52.50 
nS W530) 08) 75 0.05 52.55 
16 0.82 Uo® 105 2.85 55.40 
17 2.46 UYod 105 5.30 60.70 
18 4.10 8.5 105 3e5 64.25 
19 5.74 8.5 105 1.85 66.10 
20 7.38 Vod 105 1.20 67.30 
21 9.02 7.5 105 0.75 68.05 
22 10.66 8.5 105 0.45 68.50 
23 12.30 9.5 105 0.20 68.70 
24 13.94 10.5 105 0.10 68.80 
25 15.58 10.5 105 0.05 68.85 
26 18.04 12.0 105 0.05 68.90 
27 0.82 Uo) 135 11.05 YDoOD 
28 2.46 Tod 135 8.50 88.45 
29 4.10 Vod) 135 4.90 93.35 
30 5.74 Tod 135 2.95 96.30 
31 7.38 Vos 135 2.20 98.50 
32 9.02 Uod ILS) 0.95 99.45 
33 10.66 9.5 1335) 0.30 MoD 
34 12.30 9.5 135 0.10 99.85 
35 13.94 10.5 135 0.05 99.90 
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Figure 40. Location of entrance channel 


each cell face. The sum of the two discharges (wave-induced and tidal) 


constituted the total water discharge. 


112. The dimensionless sediment concentration for the surf zone and be- 


yond the surf zone was calculated from equations described earlier. In the 


numerical model, the sediment concentration was defined at the center of a 


cell. The water discharge was then multiplied with the sediment concentration 


to obtain sediment discharge on a cell face. 


113. For a given tidal and wave-induced flow condition the continuity 


equation (Equation 97) was solved to determine the change in the bed elevation 


Ah of a cell. A time-step At of 1 hr was used for the sediment transport 


A sensitivity analysis on time-step (At = 


model. = 10, 30, and 60 min) showed 
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that the desired accuracy could be achieved with a 1l-hr time-step. A porosity 
p of 40 percent was used for sand. 


114, The new water depth d_, of a cell was calculated by adding Ah 


WwW 
to water depth d. of the previous time-step. If Ah was negative, the cell 
eroded; and if positive, the cell accreted. The next time-step began with the 
calculation of water depth d, in each cell of the entire shore processes 
grid. This was accomplished by adding the change in tidal elevation to dj. 
d, was then used to calculate the dimensionless sediment concentrations be- 
yond the surf zone using the modified Ackers and White formulation. The surf 
zone concentrations were kept constant because the wave field was constant. 
This procedure was followed for three time-steps. 

115. Every fourth time-step a new set of records consisting of wave 
fields, wave-induced currents and setups, and tidal currents and elevations 
was read in accordance with the file sequence presented earlier. The tech- 
niques described in the previous paragraph were repeated between every fourth 
time-step. 

116. The bathymetric changes in the Oregon Inlet area were calculated for 
1 year. The location of the Oregon Inlet channel (defined earlier) at the 
beginning and end of each month was saved for analysis. In addition, the 
water depths over the entire shore processes grid were saved each month to 
determine monthly erosion and deposition. 

117. The volume of sediment trapped in the Oregon Inlet channel was com- 
puted each month for a year. Table 3 shows the monthly deposition. A total 
of 1,055,990 cu yd of materials was trapped in the Oregon Inlet channel during 
the 12-month period. A calculation using the CERC formula given in the SPM 
(1984) revealed that 1,080,000 cu yd of materials could be transported during 
a year under an average wave condition applicable to Oregon Inlet. 

118. Table 4 shows that erosion occurred between the south jetty and the 
channel boundary on the south side. A total of 660,000 cu yd of materials was 
eroded during the year. The material was either deposited in the channel or 
transported offshore by tidal and wave-induced currents. 

119. Contour plots were made at intervals of 3 months to show erosion and 
accretion in the Oregon Inlet channel. Plates 4-8 show contours of erosion 
and accretion (in feet) at the end of the Ist, 3rd, 6th, 9th, and 12th months, 
respectively. Plates 9-13 present zones of erosion and accretion in the 


channel at the end of the lst, 3rd, 6th, 9th and 12th months, respectively. 
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Table 3 


Volume of Materials Trapped in the 
Oregon Inlet Channel 


Volume 

Month cu_yd 
1 78,055 
2 100,833 
3 87,500 
4 81,667 
3) 92,778 
6 84,166 
7 86,243 
8 89,360 
9 83,068 
10 90,165 
iLaL 87,483 
2) 94,672 


Total 1,055,990 


Table 4 


Volume of Materials Eroded between the South Jetty 
and the South Channel Boundary 


Volume Eroded 
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Month cu_yd 
3335 333)3) 

2 50,833 
3 51,944 
4 55,000 
5 DoWoSS 
6 51,667 
7/ 57,248 
8 56,556 
9 60,373 
10 56,845 
11 59,778 
12 50,590 
Total 660,000 


120. Flow cross sections taken at various locations along the Oregon 
Inlet channel (Plate 13) at the beginning and the end of the l-year numerical 
simulation showed (Figure 41) that the part of the Oregon Inlet channel 
sheltered by the south jetty could move on the average about 150 ft per year 


toward the south jetty (for Figure 41, mean sea level (MSL)is considered to be 


a datum 1.0 ft above MLW). 
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PART IV: SUMMARY AND CONCLUSIONS 


121. A system of numerical models was developed to simulate coastal and 
inlet processes near Oregon Inlet, North Carolina. The system included models 
for wave propagation, wave-induced currents and setup, sediment transport 
within and beyond the surf zone, and profile response (onshore-offshore trans-— 
port). Results from a separate study on numerical simulation of tides and 
storm surge for Oregon Inlet were utilized in the present investigation. 

122. The emphasis in development of the models was on computational 
efficiency and the ability to handle complex bathymetries and large numerical 
grids encountered in practical engineering problems. 

123. The models of the system were tested individually and validated 
against available analytical solutions as well as laboratory and field data. 
All of these tests and comparisons were successful. 

124. The computational costs for the models of the system were relatively 
modest so that long-term simulations could be performed economically with the 
system. 

125. As a test for an extreme event, the Ash Wednesday storm of March 
1962 was simulated with the profile response model. The calculated erosion 
amounts of the shore-normal profiles for Bodie Island and Pea Island (on 
either side of Orégon Inlet) were compared with measured values. There was 
good agreement, especially considering the uncertainty of the measurements. 

126. As an alternative to the stabilization of the Oregon Inlet entrance 
channel by construction of two jetties, a nonstructural solution proposed by 
DOI was evaluated using the profile response model. The DOI solution involved 
disposing the dredged material from the inlet channel in the nearshore region 
with the idea that the material would be dispersed shoreward by wave action at 
a rate sufficient to prevent dredging-induced beach erosion. On the basis of 
a feasibility study conducted by SAW, two disposal schemes involving placement 
of dredged material in depths of 11 and 17 ft were simulated using wave condi- 
tions for a typical year. Various scenarios were considered, and simulations 
were performed for 1- and 5-year durations. The results of these tests indi- 
cated that on the average only 25 percent of the material disposed in the 
nearshore migrated toward the shore in a year. This was not sufficient to 


prevent dredging-induced beach erosion. An independent analysis of the 
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disposal plan carried out by COER for SAW, using an n-line model, produced 
similar results. 

127. In order to perform an ocean bar channel dredging analysis, SAW 
needed to know the period of time that dredges of the CURRITUCK and 
ATCHAFALAYA/MERMENTAU classes could operate in the entrance channel under the 
influence of waves. To study this problem, the wave propagation model was 
run, allowing for wave-current interactions. The tidal currents were calcu- 
lated by the numerical tide model. The effects of jetty spacing on wave 
height and steepness were studied for a variety of typical wave conditions. 
Using these results, SAW determined the limiting wave heights for dredging 
operations in Oregon Inlet to be deepwater significant heights of 3.0 ft for 
the CURRITUCK and 4.0 ft for the ATCHAFALAYA/MERMENTAU classes. This was con- 
firmed by SAW on review of the daily logs of the dredge HYDE which operated at 
Oregon Inlet. 

128. The modeling system was used to study the erosion and accretion in 
the entrance channel, as well as the lateral movement of the channel in the 
presence of the south jetty alone, to simulate a construction sequence in 
which the south jetty was built before the north jetty. For this purpose, an 
average year's wave climate obtained from WESWIS was modeled. The wave pro- 
pagation, wave-induced current, and longshore sediment transport models were 
run in succession. Also, information on tidal currents and elevations for the 
mean tide condition obtained from the numerical tide model were used. The 
sediment transport model simulated an average year, using a time-step of 1 hr. 
Computer plots showing contours and zones of erosion and accretion were made 
at intervals of 3 months. 

129. The results of the l-year simulation showed that a total of 
1,055,990 cu yd of materials was trapped in the entrance channel during the 
12-month period, whereas 660,000 cu yd of materials were eroded between the 
southern boundary of the channel and the south jetty. 

130. On the basis of flow cross sections taken at various locations along 
the channel at the beginning and end of the numerical simulation, it was 
determined that the entrance channel could move on the average about 150 ft/yr 


toward the south jetty. 
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APPENDIX A: NOTATION 


Mapping constants for region p in x-direction 
Mapping constants for region q in y-direction 
Area off cell 

Drag coefficient, wave celerity 

Wave celerity, coefficient 

Wave group velocity 

Chezy coefficient 

Average sediment concentration (dimensionless) 
Local water depth 

New water depth of a cell 


Sediment diameter exceeded in size by 65 percent (by weight) of 
sediment sample 


Median diameter of sediment 

Wave energy density = ogH=/8 

Wave friction factor 

Wave friction factor based on bed roughness 
Wave friction factor with D as bed roughness 
Acceleration due to gravity 

Local still-water depth 


Depth of point dividing the lower limit of D-profile and the 
upper limit of transition area 


Position of maximum runup 

Wave height 

Maximum wave height in the spectrum 
Cell indices 

Wave number 

Wave number vector 


Al 


k k Components of wave number vector 
Length of onshore profile 
Lo; Length of offshore profile 


m,n Indices for cell center 


n+1/2 Indices for cell faces 

M, N Cell indices in FORTRAN 
n Ratio of group velocity to wave celerity, C¢/C 
N Empirical mixing coefficient of Longuet—Higgins 
P Porosity of sediment 
}2 Mixing parameter of Longuet-—Higgins 

BROnMRGS Recursion poetricients 

Qy Total sediment transport rate across surf zone 
Ty Bed roughness 


s Arbitrary variable, mass density of sediment relative to that of 
fluid 


S Total sediment transport rate per unit width 


Si ij Radiation stress 
Sxx»SxySyy Radiation stresses 
Sin»Sout Total sediment transport rates into and out of a cell 
oe Sediment transport rate at position i at time t 
t Time 

T Wave period 

Ws Period of wave when there is no current 

Uy Wave orbital velocity at the bottom 
<[eorp > Time average of the absolute value of the wave orbital velocity 


at bottom 


A2 


Velocity components in x- and y-directions 
Mean velocity of flow, longshore current velocity at 
Shear velocity 

Current velocity vector 

Value of Lo; - Ly; for equilibrium 

Value of W; at still-water line 
Coordinates in real space 

Width of surf zone 

Dimensionless coordinate, x/Xp 
Dimensionless grain diameter 

Arbitrary variable 

Coordinates in computational space 

Angle beach makes with the horizontal 
Breaking index 

Difference operator, total depth of D-profile 
Change in bed elevation of a cell 
Dimensionless position in D-profile 
Time-step 

Cell dimensions in real space 

Cell dimensions in computational space 
Eddy viscosities in x- and y-directions 
Mean free surface displacement 

Angle of wave propagation 

Wave direction in deep water 

Deepwater wave length 

Grid expansion coefficients 


Kinematic viscosity of fluid 


A3 


Superscripts 


Subscripts 
b 


SJoUCMtE eS oc 


Mass density of sea water 


Dimensionless radiation stresses 


Bottom friction stresses in x- and y-—directions 
Lateral shear stress due to turbulent mixing 
Radian frequency as it appears to a stationary observer 


Radian frequency when there is no current, 20/T) 


Previous time level 
Present time level 
Next time level 


Intermediate time level 


At breaking 
Partial derivative with respect to time 
Partial derivative with respect to Oy 


Partial derivative with respect to Oh 
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